####################################
# Alternative selectors 
####################################

rm(list=ls())

#sink("~/Dropbox/Crime Chile/11_replication/15_RD_figures_alt_selectors_crime.txt")

library(Hmisc)
library(ggplot2)
library(stargazer)
library(foreign)
library(rdrobust)
library(rdd)
library(readstata13)

################
# Prepare data 
################

# reada data
d = read.dta13("~/Dropbox/Crime Chile/11_replication/local_crime_data_chile_2022january.dta")   
names(d)

############################
# msesum
############################

# Homicides
summary(rdrobust(d$homicide_rate_s,d$margin,cluster = d$cluster,all=TRUE, bwselect = "msesum"))
cairo_pdf(file="~/Dropbox/Crime Chile/11_replication/figureA7a.pdf", 
          width=7, 
          height=7)
rdplot(y = d$homicide_rate_s, x = d$margin, h= 0.183, nbins = 100, subset = -0.183 <= d$margin & d$margin <= 0.183,
        binselect="esmv", kernel="triangular", col.lines = "black", col.dots = "black", p=1, y.lim = c(-1, 1), title = "                                       Robust CI: [-0.407 , 0.172]",
        y.label = "Homicide", x.label = "Margin of victory")
dev.off() 

# Rapes
summary(rdrobust(d$rape_rate_s,d$margin,cluster = d$cluster,all=TRUE, bwselect = "msesum"))
cairo_pdf(file="~/Dropbox/Crime Chile/11_replication/figureA7b.pdf", 
          width=7, 
          height=7)
(rdplot(y = d$rape_rate_s, x = d$margin,  h= 0.117, nbins = 100, subset = -0.117 <= d$margin & d$margin <= 0.117, 
        binselect="esmv", kernel="triangular", col.lines = "black", col.dots = "black", p=1, y.lim = c(-1, 1), title = "                                       Robust CI: [-0.353 , 0.097]",
        y.label = "Rape", x.label = "Margin of victory"))
dev.off()

# Assault
summary(rdrobust(d$assault_rate_s,d$margin,cluster = d$cluster,all=TRUE, bwselect = "msesum"))
cairo_pdf(file="~/Dropbox/Crime Chile/11_replication/figureA7c.pdf", 
          width=7, 
          height=7)
(rdplot(y = d$assault_rate_s, x = d$margin,  h= 0.138, nbins = 100, subset = -0.138 <= d$margin & d$margin <= 0.138, 
        binselect="esmv", kernel="triangular", col.lines = "black", col.dots = "black", p=1, y.lim = c(-1, 1), title = "                                       Robust CI: [-0.433 , 0.109]",
        y.label = "Assault", x.label = "Margin of victory"))
dev.off()

# Theft
summary(rdrobust(d$theft_rate_s,d$margin,cluster = d$cluster,all=TRUE, bwselect = "msesum"))
cairo_pdf(file="~/Dropbox/Crime Chile/11_replication/figureA8a.pdf", 
          width=7, 
          height=7)
(rdplot(y = d$theft_rate_s, x = d$margin,  h= 0.117, nbins = 100, subset = -0.117 <= d$margin & d$margin <= 0.117, 
        binselect="esmv", kernel="triangular", col.lines = "black", col.dots = "black", p=1, y.lim = c(-1, 1), title = "                                       Robust CI: [-0.833 , -0.234]",
        y.label = "Theft", x.label = "Margin of victory"))
dev.off()

# Robbery
summary(rdrobust(d$robbery_rate_s,d$margin,cluster = d$cluster,all=TRUE, bwselect = "msesum"))
cairo_pdf(file="~/Dropbox/Crime Chile/11_replication/figureA8b.pdf", 
          width=7, 
          height=7)
(rdplot(y = d$robbery_rate_s, x = d$margin,  h= 0.129, nbins = 100, subset = -0.129 <= d$margin & d$margin <= 0.129, 
        binselect="esmv", kernel="triangular", col.lines = "black", col.dots = "black", p=1, y.lim = c(-1, 1), title = "                                       Robust CI: [-0.712 , -0.041]",
        y.label = "Robbery", x.label = "Margin of victory"))
dev.off()

# Robbery by surprise
summary(rdrobust(d$robbery_surprise_rate_s,d$margin,cluster = d$cluster,all=TRUE, bwselect = "msesum"))
cairo_pdf(file="~/Dropbox/Crime Chile/11_replication/figureA8c.pdf", 
          width=7, 
          height=7)
(rdplot(y = d$robbery_surprise_rate_s, x = d$margin,  h= 0.135, nbins = 100, subset = -0.135 <= d$margin & d$margin <= 0.135, 
        binselect="esmv", kernel="triangular", col.lines = "black", col.dots = "black", p=1, y.lim = c(-1, 1), title = "                                       Robust CI: [-0.693 , 0.009]",
        y.label = "Robbery by surprise", x.label = "Margin of victory"))
dev.off()

############################
# msetwo
############################

# Homicides
summary(rdrobust(d$homicide_rate_s,d$margin,cluster = d$cluster,all=TRUE, bwselect = "msetwo"))
cairo_pdf(file="~/Dropbox/Crime Chile/11_replication/figureA9a.pdf", 
          width=7, 
          height=7)
rdplot(y = d$homicide_rate_s, x = d$margin,  nbins = 100, subset = -0.173 <= d$margin & d$margin <= 0.151,
       binselect="esmv", kernel="triangular", col.lines = "black", col.dots = "black", p=1, y.lim = c(-1, 1), title = "                                       Robust CI: [-0.407 , 0.172]",
       y.label = "Homicide", x.label = "Margin of victory")
dev.off() 

# Rapes
summary(rdrobust(d$rape_rate_s,d$margin,cluster = d$cluster,all=TRUE, bwselect = "msetwo"))
cairo_pdf(file="~/Dropbox/Crime Chile/11_replication/figureA9b.pdf", 
          width=7, 
          height=7)
(rdplot(y = d$rape_rate_s, x = d$margin,  nbins = 100, subset = -0.202 <= d$margin & d$margin <= 0.105, 
        binselect="esmv", kernel="triangular", col.lines = "black", col.dots = "black", p=1, y.lim = c(-1, 1), title = "                                       Robust CI: [-0.285 , 0.115]",
        y.label = "Rape", x.label = "Margin of victory"))
dev.off()

# Assault
summary(rdrobust(d$assault_rate_s,d$margin,cluster = d$cluster,all=TRUE, bwselect = "msetwo"))
cairo_pdf(file="~/Dropbox/Crime Chile/11_replication/figureA9c.pdf", 
          width=7, 
          height=7)
(rdplot(y = d$assault_rate_s, x = d$margin,  nbins = 100, subset = -0.193 <= d$margin & d$margin <= 0.119, 
        binselect="esmv", kernel="triangular", col.lines = "black", col.dots = "black", p=1, y.lim = c(-1, 1), title = "                                       Robust CI: [-0.373 , 0.124]",
        y.label = "Assault", x.label = "Margin of victory"))
dev.off()

# Theft
summary(rdrobust(d$theft_rate_s,d$margin,cluster = d$cluster,all=TRUE, bwselect = "msetwo"))
cairo_pdf(file="~/Dropbox/Crime Chile/11_replication/figureA10a.pdf", 
          width=7, 
          height=7)
(rdplot(y = d$theft_rate_s, x = d$margin,  nbins = 100, subset = -0.093 <= d$margin & d$margin <= 0.152, 
        binselect="esmv", kernel="triangular", col.lines = "black", col.dots = "black", p=1, y.lim = c(-1, 1), title = "                                       Robust CI: [-0.811 , -0.163]",
        y.label = "Theft", x.label = "Margin of victory"))
dev.off()

# Robbery
summary(rdrobust(d$robbery_rate_s,d$margin,cluster = d$cluster,all=TRUE, bwselect = "msetwo"))
cairo_pdf(file="~/Dropbox/Crime Chile/11_replication/figureA10b.pdf", 
          width=7, 
          height=7)
(rdplot(y = d$robbery_rate_s, x = d$margin,  nbins = 100, subset = -0.112 <= d$margin & d$margin <= 0.241, 
        binselect="esmv", kernel="triangular", col.lines = "black", col.dots = "black", p=1, y.lim = c(-1, 1), title = "                                       Robust CI: [-0.719 , -0.026]",
        y.label = "Robbery", x.label = "Margin of victory"))
dev.off()

# Robbery by surprise
summary(rdrobust(d$robbery_surprise_rate_s,d$margin,cluster = d$cluster,all=TRUE, bwselect = "msetwo"))
cairo_pdf(file="~/Dropbox/Crime Chile/11_replication/figureA10c.pdf", 
          width=7, 
          height=7)
(rdplot(y = d$robbery_surprise_rate_s, x = d$margin,  nbins = 100, subset = -0.108 <= d$margin & d$margin <= 0.194, 
        binselect="esmv", kernel="triangular", col.lines = "black", col.dots = "black", p=1, y.lim = c(-1, 1), title = "                                       Robust CI: [-0.802 , -0.051]",
        y.label = "Robbery by surprise", x.label = "Margin of victory"))
dev.off()

##################
# TABLE A17
##################

# Homicide
summary(rdrobust(d$homicide_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msesum"))

pe1 = as.numeric(rdrobust(d$homicide_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msesum")$Estimate[1])
pv1 = as.numeric(rdrobust(d$homicide_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msesum")$pv[3])
ci1 = as.numeric(rdrobust(d$homicide_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msesum")$ci[3,])
oss1 = as.numeric(rdrobust(d$homicide_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msesum")$N[1] + rdrobust(d$homicide_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msesum")$N[2])
ess1 = as.numeric(rdrobust(d$homicide_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msesum")$N_h[1] + rdrobust(d$homicide_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msesum")$N_h[2])
bw1 = as.numeric(rdrobust(d$homicide_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msesum")$bws[1,1])

# Rape
summary(rdrobust(d$rape_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msesum"))

pe2 = as.numeric(rdrobust(d$rape_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msesum")$Estimate[1])
pv2 = as.numeric(rdrobust(d$rape_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msesum")$pv[3])
ci2 = as.numeric(rdrobust(d$rape_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msesum")$ci[3,])
oss2 = as.numeric(rdrobust(d$rape_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msesum")$N[1] + rdrobust(d$rape_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msesum")$N[2])
ess2 = as.numeric(rdrobust(d$rape_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msesum")$N_h[1] + rdrobust(d$rape_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msesum")$N_h[2])
bw2 = as.numeric(rdrobust(d$rape_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msesum")$bws[1,1])

# Assault 
summary(rdrobust(d$assault_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msesum"))

pe3 = as.numeric(rdrobust(d$assault_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msesum")$Estimate[1])
pv3 = as.numeric(rdrobust(d$assault_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msesum")$pv[3])
ci3 = as.numeric(rdrobust(d$assault_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msesum")$ci[3,])
oss3 = as.numeric(rdrobust(d$assault_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msesum")$N[1] + rdrobust(d$assault_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msesum")$N[2])
ess3 = as.numeric(rdrobust(d$assault_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msesum")$N_h[1] + rdrobust(d$assault_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msesum")$N_h[2])
bw3 = as.numeric(rdrobust(d$assault_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msesum")$bws[1,1])

# Theft
summary(rdrobust(d$theft_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msesum"))

pe4 = as.numeric(rdrobust(d$theft_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msesum")$Estimate[1])
pv4 = as.numeric(rdrobust(d$theft_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msesum")$pv[3])
ci4 = as.numeric(rdrobust(d$theft_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msesum")$ci[3,])
oss4 = as.numeric(rdrobust(d$theft_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msesum")$N[1] + rdrobust(d$theft_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msesum")$N[2])
ess4 = as.numeric(rdrobust(d$theft_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msesum")$N_h[1] + rdrobust(d$theft_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msesum")$N_h[2])
bw4 = as.numeric(rdrobust(d$theft_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msesum")$bws[1,1])

# Robbery
summary(rdrobust(d$robbery_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msesum"))

pe5 = as.numeric(rdrobust(d$robbery_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msesum")$Estimate[1])
pv5 = as.numeric(rdrobust(d$robbery_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msesum")$pv[3])
ci5 = as.numeric(rdrobust(d$robbery_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msesum")$ci[3,])
oss5 = as.numeric(rdrobust(d$robbery_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msesum")$N[1] + rdrobust(d$robbery_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msesum")$N[2])
ess5 = as.numeric(rdrobust(d$robbery_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msesum")$N_h[1] + rdrobust(d$robbery_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msesum")$N_h[2])
bw5 = as.numeric(rdrobust(d$robbery_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msesum")$bws[1,1])

# Robbery by surprise
summary(rdrobust(d$robbery_surprise_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msesum"))

pe6 = as.numeric(rdrobust(d$robbery_surprise_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msesum")$Estimate[1])
pv6 = as.numeric(rdrobust(d$robbery_surprise_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msesum")$pv[3])
ci6 = as.numeric(rdrobust(d$robbery_surprise_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msesum")$ci[3,])
oss6 = as.numeric(rdrobust(d$robbery_surprise_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msesum")$N[1] + rdrobust(d$robbery_surprise_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msesum")$N[2])
ess6 = as.numeric(rdrobust(d$robbery_surprise_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msesum")$N_h[1] + rdrobust(d$robbery_surprise_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msesum")$N_h[2])
bw6 = as.numeric(rdrobust(d$robbery_surprise_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msesum")$bws[1,1])

pe = c(pe1,pe2,pe3,pe4,pe5,pe6)
pv = c(pv1,pv2,pv3,pv4,pv5,pv6)
ci_lower = c(ci1[1],ci2[1],ci3[1],ci4[1],ci5[1],ci6[1])
ci_upper = c(ci1[2],ci2[2],ci3[2],ci4[2],ci5[2],ci6[2])
oss = c(oss1,oss2,oss3,oss4,oss5,oss6)
ess = c(ess1,ess2,ess3,ess4,ess5,ess6)
bw = c(bw1,bw2,bw3,bw4,bw5,bw6)
variable = c("Homicide","Rape","Assault","Theft","Robbery","Robbery by surprise")
data = data.frame(variable,pe,pv,ci_lower,ci_upper,oss,ess,bw)

colnames(data) = c("","Point Estimate","Robust P-value","Robust 95% Confidence Interval lower bound","Robust 95% Confidence Interval upper bound","Overall sample size","Effective sample size","MSE bandwidth")
data

stargazer(data, summary=FALSE, rownames=FALSE, out="~/Dropbox/Crime Chile/11_replication/tableA17.html")

##################
# TABLE A18
##################

# Homicide
summary(rdrobust(d$homicide_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msetwo"))

pe1 = as.numeric(rdrobust(d$homicide_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msetwo")$Estimate[1])
pv1 = as.numeric(rdrobust(d$homicide_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msetwo")$pv[3])
ci1 = as.numeric(rdrobust(d$homicide_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msetwo")$ci[3,])
oss1 = as.numeric(rdrobust(d$homicide_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msetwo")$N[1] + rdrobust(d$homicide_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msetwo")$N[2])
ess1 = as.numeric(rdrobust(d$homicide_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msetwo")$N_h[1] + rdrobust(d$homicide_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msetwo")$N_h[2])
bw1 = as.numeric(rdrobust(d$homicide_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msetwo")$bws[1,1])

# Rape
summary(rdrobust(d$rape_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msetwo"))

pe2 = as.numeric(rdrobust(d$rape_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msetwo")$Estimate[1])
pv2 = as.numeric(rdrobust(d$rape_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msetwo")$pv[3])
ci2 = as.numeric(rdrobust(d$rape_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msetwo")$ci[3,])
oss2 = as.numeric(rdrobust(d$rape_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msetwo")$N[1] + rdrobust(d$rape_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msetwo")$N[2])
ess2 = as.numeric(rdrobust(d$rape_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msetwo")$N_h[1] + rdrobust(d$rape_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msetwo")$N_h[2])
bw2 = as.numeric(rdrobust(d$rape_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msetwo")$bws[1,1])

# Assault 
summary(rdrobust(d$assault_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msetwo"))

pe3 = as.numeric(rdrobust(d$assault_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msetwo")$Estimate[1])
pv3 = as.numeric(rdrobust(d$assault_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msetwo")$pv[3])
ci3 = as.numeric(rdrobust(d$assault_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msetwo")$ci[3,])
oss3 = as.numeric(rdrobust(d$assault_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msetwo")$N[1] + rdrobust(d$assault_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msetwo")$N[2])
ess3 = as.numeric(rdrobust(d$assault_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msetwo")$N_h[1] + rdrobust(d$assault_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msetwo")$N_h[2])
bw3 = as.numeric(rdrobust(d$assault_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msetwo")$bws[1,1])

# Theft
summary(rdrobust(d$theft_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msetwo"))

pe4 = as.numeric(rdrobust(d$theft_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msetwo")$Estimate[1])
pv4 = as.numeric(rdrobust(d$theft_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msetwo")$pv[3])
ci4 = as.numeric(rdrobust(d$theft_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msetwo")$ci[3,])
oss4 = as.numeric(rdrobust(d$theft_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msetwo")$N[1] + rdrobust(d$theft_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msetwo")$N[2])
ess4 = as.numeric(rdrobust(d$theft_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msetwo")$N_h[1] + rdrobust(d$theft_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msetwo")$N_h[2])
bw4 = as.numeric(rdrobust(d$theft_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msetwo")$bws[1,1])

# Robbery
summary(rdrobust(d$robbery_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msetwo"))

pe5 = as.numeric(rdrobust(d$robbery_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msetwo")$Estimate[1])
pv5 = as.numeric(rdrobust(d$robbery_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msetwo")$pv[3])
ci5 = as.numeric(rdrobust(d$robbery_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msetwo")$ci[3,])
oss5 = as.numeric(rdrobust(d$robbery_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msetwo")$N[1] + rdrobust(d$robbery_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msetwo")$N[2])
ess5 = as.numeric(rdrobust(d$robbery_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msetwo")$N_h[1] + rdrobust(d$robbery_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msetwo")$N_h[2])
bw5 = as.numeric(rdrobust(d$robbery_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msetwo")$bws[1,1])

# Robbery by surprise
summary(rdrobust(d$robbery_surprise_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msetwo"))

pe6 = as.numeric(rdrobust(d$robbery_surprise_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msetwo")$Estimate[1])
pv6 = as.numeric(rdrobust(d$robbery_surprise_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msetwo")$pv[3])
ci6 = as.numeric(rdrobust(d$robbery_surprise_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msetwo")$ci[3,])
oss6 = as.numeric(rdrobust(d$robbery_surprise_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msetwo")$N[1] + rdrobust(d$robbery_surprise_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msetwo")$N[2])
ess6 = as.numeric(rdrobust(d$robbery_surprise_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msetwo")$N_h[1] + rdrobust(d$robbery_surprise_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msetwo")$N_h[2])
bw6 = as.numeric(rdrobust(d$robbery_surprise_rate_s,d$margin,all=TRUE,cluster = d$cluster,  bwselect = "msetwo")$bws[1,1])

pe = c(pe1,pe2,pe3,pe4,pe5,pe6)
pv = c(pv1,pv2,pv3,pv4,pv5,pv6)
ci_lower = c(ci1[1],ci2[1],ci3[1],ci4[1],ci5[1],ci6[1])
ci_upper = c(ci1[2],ci2[2],ci3[2],ci4[2],ci5[2],ci6[2])
oss = c(oss1,oss2,oss3,oss4,oss5,oss6)
ess = c(ess1,ess2,ess3,ess4,ess5,ess6)
bw = c(bw1,bw2,bw3,bw4,bw5,bw6)
variable = c("Homicide","Rape","Assault","Theft","Robbery","Robbery by surprise")
data = data.frame(variable,pe,pv,ci_lower,ci_upper,oss,ess,bw)

colnames(data) = c("","Point Estimate","Robust P-value","Robust 95% Confidence Interval lower bound","Robust 95% Confidence Interval upper bound","Overall sample size","Effective sample size","MSE bandwidth")
data

stargazer(data, summary=FALSE, rownames=FALSE, out="~/Dropbox/Crime Chile/11_replication/tableA18.html")

sink()
